%% 加载海洋信道数据
clear;
close all;
%% channel 1 
% load("raychan.mat");
% At = 100;
% raychan=raychan*At;
% % mesh(abs(raychan));

%% --------------------------------------------------------------------------
%
% This code forms a simulation of a wideband wireless communications system
% in multipath fading. It simulates: OFDM,  OTFS, coded-OFDM and coded-OTFS
% 
% The following files are required:
% dataGen.m
% multipathChannel.m
% modOFDM.m
% demodOFDM.m
% ISFFT.m
% SFFT.m
% equaliser.m
% plotGraphs.m
%%
%--------------------------------------------------------------------------
%
% Author: Bradley Bates
% University of Bristol, UK
% email address: bb16177@bristol.ac.uk
% May 2020
%
% Copyright (c) 2020, Bradley Bates
%
%--------------------------------------------------------------------------

%% --------------------------------------------------------------------------
%% 定义信道参数
%--------------------------------------------------------------------------
global  lambda B c
M_mod = 4;                         % 16QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
Mp = 0;
% 
% alphabet = qammod((0:2^k-1),2^k,'gray');            % alphabet of M-QAM modulation
% P_sym = mean(abs(alphabet).^2);                               % average power of alphabet for power normalization

% cpSize = 0.05;                  % 循环前缀比例

B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 

codeRate = 1/2;                 % LDPC码率
maxIterations = 20;             % 最大迭代次数LDPC decoder
totalBits = 1e6;                % The approx. total no. of bits simulated

repeats = 1;                    % 重复次数


fc = 12500;
c = 1500;                                                     % speed of light (m/s)
lambda = c/fc; 



eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
SNR_dB = 5:5:25;
SNR = 10.^(SNR_dB/10);
noise_var_sqrt = sqrt(1./SNR);
sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;                   % 噪声功率

%--------------------------------------------------------------------------
%%                    Initialise Simulation Components
%--------------------------------------------------------------------------

% 初始化OFDM参数
numSC = M; % 计算载波数量
% cpLen = floor(cpSize * numSC);    % 计算循环前缀长度
numDC = (numSC - Mzp-Mp);             % 计算数据长度

% 初始化高斯信道和误码率存储（设置为0）对象
awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');
errorRate = comm.ErrorRate('ResetInputPort',true);
errorRate1 = comm.ErrorRate('ResetInputPort',true);
errorRate2 = comm.ErrorRate('ResetInputPort',true);

% 初始化LDPC对象
parityCheck_matrix = dvbs2ldpc(codeRate);               % 创建 DVB-S.2 校验矩阵
ldpcEncoder = comm.LDPCEncoder(parityCheck_matrix);     % 编码对象
ldpcDecoder = comm.LDPCDecoder(parityCheck_matrix);     % 解码对象
ldpcDecoder.MaximumIterationCount = maxIterations;      % 分配最大迭代次数
noCodedbits = size(parityCheck_matrix,2);               % 码字比特长度

% 创建误码率储存对象
berOFDM = zeros(length(SNR),3); 
berCOFDM = zeros(length(SNR),3);
berOTFS = zeros(length(SNR),3);
berCOTFS_MRC = zeros(length(SNR),3);
berCOTFS_LMMSE = zeros(length(SNR),3);
berchannal = zeros(length(SNR),3);
errorStats_coded1 = zeros(1,3);
errorStats_coded2 = zeros(1,3);
errorStats_coded3 = zeros(1,3);
errorStats_coded4 = zeros(1,3);
errorStats_coded5 = zeros(1,3);
errorStats_coded6 = zeros(1,3);
%% channel 2
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fn2 = dftmtx(M);       % Generate the DFT matrix
Fn2 = Fn2./norm(Fn2);    % normalize the DFT matrix

Md = M-Mzp-Mp;

data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
delta_f = scs;
T = 1/delta_f;
% length_CP=cpLen;

max_speed=20;  % km/hr
[chan_coef,delay_taps,Doppler_taps,taps]=Gen_delay_Doppler_channel_parameters(N,M,fc,delta_f,T,max_speed,1);
L_set=unique(delay_taps);

[G,gs]=Gen_time_domain_channel(N,M,taps,delay_taps,Doppler_taps,chan_coef);



%% Generate delay-time and delay-Doppler channel vectors from the time domain channel.
%         [nu_ml_tilda]=Gen_delay_time_channel_vectors(N,M,l_max,gs);
[nu_ml_tilda,nu_ml,K_ml]=Gen_DT_and_DD_channel_vectors(N,M,L_set,gs); 
%% Generate block-wise time-frequency domain channel
[H_tf]=Generate_time_frequency_channel_ZP(N,M,gs,L_set);
% l_max=max(L_set);
        
%%channel output                  
% gs=Gen_time_domain_channel_OTFSvariants(N,M,delay_taps,Doppler_taps,chan_coef,length_CP,'ZP');
% r=Gen_channel_output_OTFSvariants(N,M,L_set,gs,s,sigma_2(iesn0),length_CP,variant);              
      
 %% channel 3
% tau = [0.0004	0.0024  0.0067	0.0122	0.0223	0.0348]; % Path delays in seconds
% pdb = [-13.582 -23.040 -24.505 -29.816	-38.278	-46.364]; % Path power in dB
% fs = scs;
% fd = 1.0755*car_fre/c;
% % Create a Rayleigh fading channel object
% rayleighchan = comm.RayleighChannel('SampleRate', fs, 'PathDelays', tau, 'AveragePathGains', pdb, 'MaximumDopplerShift', fd);
% % % Pass the upsampled signal through the channel
% % fadedSignal = rayleighchan(upsampledSignal);


%% 迭代
for repetition=1:repeats                       % 迭代次数
    
    %生成数据
    [dataIn, dataBits_in, codedData_in, packetSize, numPackets, numCB] = dataGen(k,numDC,ofdmSym,totalBits,codeRate,ldpcEncoder);
     data_na = reshape(dataIn',k*numDC*ofdmSym,packetSize).';
    
%     % Generate Rayleigh Fading Channel Impulse Response
%     txSig_size = zeros((numSC+cpLen),ofdmSym); 
%     Assign the size of the channel
%     rayChan = multipathChannel(cpSize, scs, txSig_size, velocity);   % Create fading channel impulse response   
%%                       OTFS BER Calculation
%--------------------------------------------------------------------------
    
    % 计算信噪比
  


    % 不同EbNo循环
    for m = 1:length(SNR)
        % 不同包循环
        for j = 1:numPackets
            rxframeBuffer1 = [];                 % Initialise matrix 
            rxframeBuffer2 = [];  
            % 不同帧循环
            for u = 1:packetSize
                %一次发送一个帧\
   
               
             data=qammod(reshape(data_na(u,:)',k,Md*N), M_mod,'gray','InputType','bit');
             X = Generate_2D_data_grid(N,M,data,data_grid);
        
        
        %% OTFS modulation%%%%
            X1_tilda=X*Fn';
            X2_tilda=Fn2'*X;
            s1= reshape(X1_tilda,N*M,1);
            s2= reshape(X2_tilda,N*M,1);
          
        
         %% channel output%%%%%

        r1=zeros(N*M,1);
        r2=zeros(N*M,1);
        noise1= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));
        noise2= sqrt(sigma_2(m)/2)*(randn(size(s2)) + 1i*randn(size(s2)));
        l_max=max(delay_taps);
        for q=1:N*M
            for l=(L_set+1)
                if(q>=l)
                    r1(q)=r1(q)+gs(l,q)*s1(q-l+1);
                    r2(q)=r2(q)+gs(l,q)*s2(q-l+1);
                end
            end
        end
        r1=r1+noise1;
        r2=r2+noise2;
         %% OTFS demodulation%%%%
        Y1_tilda=reshape(r1,M,N);
        Y1 = Y1_tilda*Fn;
        
        Y2_tilda=reshape(r2,M,N);
        Y2 = Fn2*Y2_tilda;
        y1=reshape(Y1.',N*M,1);
        y2=reshape(Y2.',N*M,1);
            %% detector compare%%%%
%             n_ite_MPA = 15;     % maximum number of MPA detector iterations
                n_ite_MRC = 40;     % maximum number of MRC detector iterations 
            % damping parameter - reducing omega improves error performance at the cost of increased detector iterations
                omega = 1;
            if M_mod >= 64
                omega = 0.25;   % set omega to a smaller value (for example: 0.05) for modulation orders greater than 64-QAM
            end
            decision = 1;       % 1-hard decision, 0-soft decision
            init_estimate = 1;  % 1-use the TF single tap estimate as the initial estimate for MRC detection, 0-initialize the symbol estimates to 0 at the start of MRC iteration

            [est_info_bits_MRC_es1,~,~] = MRC_delay_time_detector(N,M,Md,M_mod,sigma_2(m),data_grid,r1,H_tf,nu_ml_tilda,L_set,omega,decision,init_estimate,n_ite_MRC); 
            [est_info_bits_MRC_es2,~,~] = MRC_delay_time_detector(N,M,Md,M_mod,sigma_2(m),data_grid,r2,H_tf,nu_ml_tilda,L_set,omega,decision,init_estimate,n_ite_MRC); 
%             [est_info_bits_LMMSE,~,~]= Algorithm1_detector(N,M,Md,M_mod,sigma_2(m),data_grid,rb,H_tf_es,L_set,omega,decision,init_estimate,n_ite_MRC,K_ml,rxSubframe');
            rxframeBuffer1 = [rxframeBuffer1;est_info_bits_MRC_es1];     % 储存帧信息
            rxframeBuffer2 = [rxframeBuffer2;est_info_bits_MRC_es2];  
            end

            dataOut1 = rxframeBuffer1;
            dataOut2 = rxframeBuffer2;
%             errorStats_channal = errorRate2(dataIn,dataOut1,0);
            % 一个包未编码信息计算误码率
%           dataOut = qamdemod(qamRx,M,'OutputType','bit','UnitAveragePower',true);% QAM解调
%             codedData_out1 = randdeintrlv(dataOut1,4831);                            % 去交织
%             codedData_out1(numel(codedData_in)+1:end) = [];                         % 删除填充比特
%             errorStats_uncoded = errorRate(codedData_in,codedData_out1,0);          % 收集错误信息
            

            % 一个包编码信息计算误码率
            
            codedData_out1 = randdeintrlv(dataOut1,4831);                      %去交织
            codedData_out1(numel(codedData_in)+1:end) = [];                   % 删除填充比特
            
            codedData_out2 = randdeintrlv(dataOut2,4831);                      %去交织
            codedData_out2(numel(codedData_in)+1:end) = [];                   % 删除填充比特
            
            % 解码
            dataBits_out1 = [];  
            dataBits_out2 = []; 
%             cc=ldpcEncoder(dataBits_in);
%             dd=ldpcDecoder(cc);
            
            dataOut_buffer1 = codedData_out1;
            dataOut_buffer2 = codedData_out2;
            for q = 1:numCB
                dataBits_out1 = [dataBits_out1;ldpcDecoder(1-2*dataOut_buffer1(1:noCodedbits))]; %解码
                dataOut_buffer1(1:noCodedbits) = [];     
                
                dataBits_out2 = [dataBits_out2;ldpcDecoder(1-2*dataOut_buffer2(1:noCodedbits))]; %解码
                dataOut_buffer2(1:noCodedbits) = [];  
            end
            errorStats_coded1 = errorRate1(dataBits_in,dataBits_out1,0);     % 计算错误信息
            errorStats_coded2 = errorRate2(dataBits_in,dataBits_out2,0);
            
        end
        berCOTFS_MRC(m,:) =berCOTFS_MRC(m,:)+ errorStats_coded1.';                                  % 保存未编码错误信息
        berCOTFS_LMMSE(m,:) = berCOTFS_LMMSE(m,:)+errorStats_coded2.';  
        % 保存编码错误信息

        errorStats_coded1 = errorRate1(dataBits_in,dataBits_out1,1);          % 重置
        errorStats_coded2 = errorRate2(dataBits_in,dataBits_out2,1);
    end
    
end
%% 计算误码率
for i = 1:size(SNR,2)
    berCOTFS_MRC(i,1) = berCOTFS_MRC(i,2)/berCOTFS_MRC(i,3);
    berCOTFS_LMMSE(i,1) = berCOTFS_LMMSE(i,2)/berCOTFS_LMMSE(i,3);
end
%--------------------------------------------------------------------------
%%                           Figures
%-------------------------------------------------------------------------- 

% 画图
% plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, SNR);
plotGraphs1(berCOTFS_MRC, berCOTFS_LMMSE, M_mod ,SNR_dB);
% save('data.mat')
